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1. Introduction 



The theory of the strong interactions, also called Quantum Chromodynamics (QCD), describes the inter- 
actions between quarks and gluons and is responsible for the existence of hadrons. Lattice-regularized 
QCD allows for the description of low-energy properties and other nonperturbative phenomena in QCD 
and has the salient property that it can be systematically improved towards the continuum limit. In lat- 
tice QCD, space-time is discretized and the functional integral of the quantum field theory is performed 
by a Markov- chain Monte-Carlo method. 

An important subject of study is the behavior of QCD in an environment exhibiting an abundance of 
particles over anti-particles. Such conditions arise, e.g., in astrophysical objects like neutron stars, and 
can be reproduced in heavy-ion collision experiments. Part of the interest comes from the existence of 
various phases in QCD, which are usually exemplified by means of the QCD phase diagram [1]. To study 
QCD at nonzero baryon density a quark chemical potential is introduced in the QCD Lagrangian. (In the 
following, we will omit the qualifier "quark" and only speak of a chemical potential.) In the presence of a 
chemical potential the QCD Dirac operator is no longer anti-Hermitian, i.e., its eigenvalues spread into 
the complex plane and its determinant will generically be complex. 

In lattice QCD the effect of dynamical fermions can be integrated out, leaving behind the determi- 
nant of the Dirac operator. Dynamical lattice simulations for QCD at nonzero chemical potential are 
problematic because the fermion determinant is complex and hence its real part can be negative, which 
prohibits its incorporation in the weight of the Monte-Carlo sampling. This is the so-called sign problem, 
which also occurs in other theories and has been the subject of a large number of investigations in recent 
years (for an incomplete list see, e.g., Refs. [2, 3, 4, 5, 6, 7, 8, 9]). While many of these works are concerned 
with a solution of the sign problem by various clever ideas, here we concentrate on an analytical study of 
the sign problem, in the hope that the results we derive will contribute to its solution. The severeness of 
the sign problem depends on the magnitude of the chemical potential, and it is therefore illuminating to 
investigate the relation between the phase factor of the determinant and the chemical potential. 

Chiral random matrix theory (chRMT) is a useful auxiliary in the study of the spectral properties of 
the Dirac operator in QCD [10, 11, 12]. Indeed, to leading order in the e-regime of QCD the spectral prop- 
erties of the Dirac operator are universal and can be described by chRMT [13] . In the presence of a chem- 
ical potential this correspondence is still valid even though the Dirac operator is now non-Hermitian. 
Appropriate random matrix models have been developed [14, 15, 16, 17], and their correspondence with 
QCD at nonzero chemical potential was verified successfully, see Ref. [18] for a review. The agreement 
of the microscopic spectral properties of the Dirac operator with the predictions of chiral random ma- 
trix theory has been confirmed for quenched lattice QCD simulations with chemical potential using the 
staggered operator [19], and more recently using the overlap operator [20, 21]. The latter operator has the 
interesting property that it satisfies the Ginsparg-Wilson relation and the trace anomaly at finite lattice 
spacing and can therefore have exact zero modes [22, 23, 24, 25, 26, 27] . This allowed us to verify the pre- 
dictions of chRMT at nonzero chemical potential for both zero and nonzero topology. The comparison 
between lattice QCD and chRMT also allows for a determination of the low-energy constants I and F of 
chiral perturbation theory. 

Motivated by this agreement, we expect that a study of the behavior of the fermion determinant in 
chRMT will give us, in certain well-defined limits, important information about the sign problem that 
is encountered in dynamical QCD simulations at nonzero chemical potential. In Ref. [28] Splittorff and 
Verbaarschot derived a solution for the average phase factor of the determinant in the microscopic limit 
of QCD (see Sec. 3.2 for a description of this limit) for the case of trivial topology using chRMT at nonzero 
chemical potential. However, to compare the overlap data of Ref. [20] with chRMT one also needs the 
RMT predictions for the average phase factor for general topology. The derivation of a formula for general 
topology is the main goal of this paper. As will be seen, the final expression contains two distinct parts. 
The first part is the generalization of the integrals representing the solution in Ref. [28] from zero to 
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arbitrary topology. The second part is a low-degree bivariate polynomial in mass and chemical potential 
which is absent for topological charge v = 0. For v ^ it gives an important contribution to the average 
phase factor of the fermion determinant, especially for small mass. As the mass goes to zero only this 
term remains and completely determines the value of the average phase. 

An important ingredient of the derivation is the ability to write the phase factor of the determinant 
as a ratio of characteristic polynomials. This quantity is recurrent in random matrix studies, both for 
theories with real [29, 30] and with complex eigenvalues [31], and its average can be computed in terms 
of Cauchy transforms of the orthogonal polynomials of the theory. To determine the phase factor of the 
determinant, the relevant Cauchy transform was computed in Ref. [28] and expressed in terms of one- 
dimensional integrals for zero topology, i.e., for square random matrices. In the present paper we extend 
the solution of the Cauchy transform to the case of rectangular matrices. This solution could also be 
relevant for other applications, like those involving time series, where one matrix dimension is typically 
much larger than the other. 

The structure of this paper is as follows. In Sec. 2 we describe the chiral random matrix model at 
nonzero chemical potential. In Sec. 3 we show how the microscopic limit of the phase of the fermion 
determinant, in both the quenched and the unquenched case, can be formally computed for such a 
matrix model in terms of a complex Cauchy transform integral. This two-dimensional integral is strongly 
oscillating, and in Sec. 4 we apply and extend the method of Ref. [28] to transform this integral into 
a much simpler and better behaved expression, involving only one-dimensional integrals and a short 
double sum (or bivariate polynomial). Explicit results for the quenched and unquenched cases as well as 
for the chiral and thermodynamic limits are given in Sec. 5. In that section we also verify the analytical 
predictions for the quenched case by random matrix simulations for various values of the topological 
charge. We conclude in Sec. 6. A number of technical details are worked out in several appendices. 



2. Non-Hermitian chiral random matrix model 



To leading order in the e-regime of QCD the spectral properties of the Dirac operator can be described by 
chRMT. In the presence of a chemical potential fi the Dirac operator D is no longer anti-Hermitian, and 
in the non-Hermitian chiral random matrix model introduced by Osborn [17] it takes the form 



DQi) = 



+ 
j'<P t + pT' t 



(2.1) 



where the matrices <J> and T 1 are complex random matrices of dimension {N+ v) x N, distributed accord- 
ing to the Gaussian weight function 



w(X) = (N/n) mN+v) exp(-JVtrX + X) . 



(2.2) 



For a detailed analysis of this model, see also Ref. [32] . For the conversion of random matrix units to 
physical units, see the beginning of Sec. 3.2. 

The parameter N will be taken to infinity when computing the microscopic limit (see Sec. 3.2). The 
matrix in Eq. (2.1) has |v| exact zero modes, which allows us to identify v with the topological charge. In 
the following, we keep v fixed as N — * oo and assume without loss of generality that v > 0. (For v < we 
can simply replace v by | v\ in the analytical results that will be computed below in the large- N limit.) The 
nonzero eigenvalues of D[^i) come in N pairs (zfc, -Zk). F° r M - 0, the are purely imaginary. 

For fixed v, the partition function of the random matrix model is given by 

r N f 
Zv> r (ju;{m/})= I dQHFV w{®)w?i>) J] det{D[fi) + m f ), (2.3) 
J f=i 
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where the integration measure is defined by 

N+v N 

dX = II 11 dReX ke dlmX k e, (2.4) 
k=i e=i 

Nf is the number of dynamical quarks, and the m f are the quark masses. The quenched case corre- 
sponds to Nf = 0, i.e., the fermion determinants are absent. 

To perform the integration over <J> and it is convenient to go to an eigenvalue representation of the 
random matrix D (p) . As shown in Ref. [17], the partition function can be rewritten, up to a normalization 
constant that depends on /i and v, as an integral over the z k , 



N N f 

/= 

where we introduced a = p 2 , the integrals over the z k are over the entire complex plane, 



n r N Nf 

Z v f (a;{m f })= H d 2 z k w v {z k ,z* k ;a)\A N {{z 2 })\ 2 l\{m 2 - z 2 k ) , (2.5) 



A N az 2 })=l\iz 2 -z 2 ) (2.6) 

k>e 

is a Vandermonde determinant, the weight function is given by 

w {z,z ;a) = \z\ exp — — (z +z )\K V \ — — — \z\ I, (2.7) 

and K v is a modified Bessel function. The quenched partition function will be denoted by Z v [a). 
The ensemble average of an observable © is given by 



I r N N f 

«®>v,N f = ^r Ud 2 z k w v {z k ,z* k ;a)\A N ({z 2 })\ 2 Yl{m 2 -z 2 k )&(z 1 ,...,z N ). (2.8) 

When there is no danger of confusion we will omit one or both of the subscripts on (6) . 

Our derivation will follow the general line of arguments given in Ref. [28] for v = 0, with the necessary 
generalizations to arbitrary topology v. To analyze the spectral properties of the random matrix model it 
is useful to introduce the orthogonal polynomials corresponding to the weight function (2.7) [17], 

l-a\ k . „( Nz 2 



where iXXz) is the generalized or associated Laguerre polynomial of order v and degree k. These orthog- 
onal polynomials satisfy the orthogonality relation 

/ d 2 zw v {z,z*\a)p v k {z;a)p v Az;a)* = rl(a)8 ke (2.10) 
Jc 

with norm 

7ra(l + a) 2fc+v A:!(fc + v)! 
W= ]^-2 ■ 

For later use we also introduce the Cauchy transform of the orthogonal polynomials defined by 

f d 2 z 

h v dm;a) = I —„ =w v {z,z* ;a)pl{z;a)* . (2.12) 

K Jc z z -m z K 
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3. Phase factor of the fermion determinant 



3. 1 The phase factor as a complex Cauchy transform 

The Dirac operator describing a massive fermion is defined as D{m; /i) = D{p) + ml, where we assume 
that m is real. If we write its determinant as detD(m;/i) = re , the phase factor can be extracted by 
forming the ratio 

e2W = det(D( M ) + m) = " m 2 -z 2 k 
det(D+ (ju) + m) l = \ m 2 - z* 2 ' 

In this expression, m is the mass of a valence quark. From the physics point of view, an interesting quan- 
tity is the ensemble average of e 2,e with two light dynamical quarks that have the same mass as the 
valence quark. This quantity tells us how the two-flavor determinant in the weight function oscillates. 
For simplicity, we shall refer to e 2 ' 8 as the phase factor of the determinant, even though it is really the 
phase of the square of the determinant. 

Because of the symmetries of (2.2), each matrix appears in the ensemble average with the same 
probability as its Hermitian conjugate. As the corresponding determinants are complex conjugate, the 
ensemble average of the phase factor is real. For strongly oscillating determinants the average phase 
factor will be close to zero, and the sign problem will be severe. On the other hand, for values of the 
chemical potential for which the average phase factor is close to unity one should still be able to perform 
dynamical simulations. 

For each topic that is treated here and in the following sections, we will first address the quenched 
case and then generalize to the unquenched case. The virtue of this approach is that the quenched case 
already contains the essential ingredients, but the arguments and the notation can be kept simple. The 
generalization to the unquenched case is straightforward but leads to somewhat more complicated ex- 
pressions. 

3.1.1 Quenched case 

The quenched ensemble average for the phase factor is given by 

2iB\ I det(Z)(jii) + m) \ Zy {a,m) 



x /N f=° \det(Dt( M ) + m)/ Wf=0 ZAa) 



i r N 

= — n 

Z v {a) Jc£=i 



N f =0 

2 2 
171 — Z 

d 2 z k w v {z k ,z* k ;a) \A N ({z 2 })\ 2 — £ . (3.2) 

m 2 - z* 2 

k 



The quantity (a, m) is the partition function of a random matrix model with one fermionic quark 
and one conjugate bosonic quark, see Ref. [28] for a detailed discussion. 

Using the formalism developed in Refs. [33, 31], the quenched average of ratios of characteristic 
polynomials can be written in terms of the orthogonal polynomials (2.9) and their Cauchy transforms 
(2.12). Applying this formalism to the quenched average phase factor (3.2) gives the compact expression 



(e 2ie \ = - 



h v N _ x [m;a) h v N {m;a) 
p^_ x (m;a) p v N {m;a) 



(3.3) 



which is a complex integral over the orthogonal polynomials due to the Cauchy transforms. This ex- 
pression (and its analog for the unquenched case, see Sec. 3.1.2) will prove to be very useful to compute 
the phase factor. Inserting the Cauchy transform (2.12) in Eq. (3.3) yields an integral over orthogonal 
polynomials, 



Pjv_i(z*;a) P v N {z*\a) 
p^_j(m;a) p v N {m\a) 



(3.4) 
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where we also used p k iz)* = p v k {z*). The well-known recurrence relation for the generalized Laguerre 
polynomials, 



[k + l)L v k+1 (*) = (fc + 1 + v)L v k {x) - xV k + 1 (x) , 
translates into a recurrence relation for the p k defined in Eq. (2.9), 



(3.5) 



1 - a 



p v k+ j (z; a) = {k + 1 + v) | \p v k {z;a) + z 2 p v k +l (z; a) . 



(3.6) 



Since the determinant remains unchanged when forming linear combinations of its columns, we can 
rewrite the phase factor (3.4) using the recurrence (3.6) as 



(e 2W ) 



d 2 z 
c z 2 - m 2 



w v (z,z*;a) 



p v N _^z*;a) z* 2 p v N + _\{z*;a) 
p^_j(m;a) m 2 p v +\{m; a) 



(3.7) 



where all the orthogonal polynomials are now of equal degree. 



3.1.2 Unquenched case 



In the presence of Nf dynamical fermion flavors with masses mi,..., rriN f , the phase factor for a valence 
quark of mass m is given by 



2W _ I det(D(p) + m) \ 
{S )N f~\det(DHp) + m)/ Nf 



i r N 

— In 

i: Jc fc=i 



m 2 -zl % 



Z v (a;{nif}) 



d 2 z k w v (z k , z* k ;a) |A w ({z 2 })| 2 — f- \\ (m 2 - z 2 ) , 



m2 - z l 2 f^"" 



(3.8) 



Nt 

where Z v [a;{m f}) is given by Eq. (2.5) . This can be written as a ratio of two partition functions, 



Nt 

Z v ' {a;{mf}) 



(3.9) 



JV f +l|l* 

where, in analogy to Eq. (3.2), Z v is the chRMT partition function with Nf + 1 fermionic quarks and 
one conjugate bosonic quark. Both partition functions can be computed using the results of Ref. [31], 

JVf+l|l* Nf 

but to apply these results we need to change the normalization and divide both Z v and Z v by 

the quenched partition function. Then the partition functions can be interpreted as averages of ratios of 
characteristic polynomials in the quenched ensemble. Applying the formalism of Ref. [31] to the numer- 
ator of Eq. (3.9) (with modified normalization) we find 



Z^ +1|1 *(a,m;{m /} ): 



r^^ANf+iim 2 ,^ 2 }) 



h v N _ x {m;a) h v N {m;a) ■■■ h v N+Nf {m;a) 
p v N _^rn;a) p v N {m;a) ••• p v N+Nf {m;a) 
P#_i(wi;a) p v N {mi;a) ••• p v N+N {mi;a) 



P N _iim Nf ;a) p v N {m Nf ;a) ••• p v N+N (m Nf ;a) 



(3.10) 
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where the matrix in the determinant is of size [Nf + 2) x (Nf + 2) . Substituting the definition of the Cauchy 
transform (2.12) gives 



Z v {a,m;{mf}) 



r d z z 

\ J z 2 -m 



■w v {z,z*;a) 



P V N+Nf {z*;a) 



r v N _ x [a) A Nf+ i (m 2 , {mj}) 

p^iz*;^ p v N {z*;a) 

P N ^ im;a) p v N {m;a) ••• p v N+Nf {m;a) 

p v N l {mi\a) p^im^a) ••• p v N+Nf (mi,a) 

P^_! ("%;«) p v N (m Nf ;a) ••• p v N+Nf {m Nf ;a) 
With the recurrence relation (3.6) for the orthogonal polynomials this becomes 

1 f d 2 z 



Z v (a,m;{mf}) = 



r] II _ l (a)A Nf+ i{m z ,{m z f }) J z z - m z 

p^tf-.a) z* 2 p v N + _\(z*;a) 
p^^imw) m 2 p v +\{m;a) 



■w v (z, z*;a) 



im 2 ) N f^p v N + _ N / + \m;a) 



v+Ne+l 



P v N _i(mi;a) mfp v ^_\{mi;a) ••• (raf) N f + p N _{ (mi;a) 



Px^irriN^a) m^^^rriN^a) ••• {m 2 Nf ) N f +1 p V ^/ +1 (m Nf ;a) 
The denominator in Eq. (3.9) (with modified normalization) can be written as [31] 



N f 

Z v {a;{rrif}) - 



1 



An f (imp 



p v N {mi;a) p v N+l {mr,a) ••• p v N+Nf _ l {m l ;a) 
p v N {m 2 ;a) p v N+l (m 2 ;a) ••• p^ +Nf _^m 2 ; a) 



p v N [m Nf ;a) p v N+l {m Nf ;a) ■■■ p v N+Nrl {m Nf ;a) 
and using the recurrence relation (3.6) for the orthogonal polynomials this can be rewritten as 



Z v f {a;{mf}) ■ 



A Nf {{m 2 }) 



p v N [m\;a) m\p v ^ ~ l (mi; a) 
p v N {m 2 ;a) m\p^ l (m 2 ;a) 



p v N {m Nf ;a) m 2 N p v + l {m Nf ;a) 



(m 2 ) N f- l p N ~ J "(mi;a) 



v+Nf-l , 



(m 2 ) N f- l p N " J '{m 2 ;a) 



V+Nr-l 



{m 2 Nf ) N f l p V * f l {m Nf ;a) 



3.2 Microscopic limit 



(3.11) 



(3.12) 



(3.13) 



(3.14) 



Universal results, i.e., results that also apply to QCD, can be obtained from chRMT in the so-called mi- 
croscopic regime. This regime is obtained by taking N — ► oo while keeping the rescaled parameters 
m - 2Nm, rhf - 2Nrrif, and a = 2Na fixed, and rescaling the spectrum using z = 2Nz. The rescaled 
random matrix parameters can be converted to the physical parameters z, m, and p using the relations 
z = zVL, m = mVE, and a - p 2 - p 2 F 2 V, where V is the four-volume. 1 Furthermore, the pion mass 
m n can be introduced through the combination p 2 lm 2 - p? 12m, where we have used the Gell-Mann- 
Oakes-Renner relation m 2 F 2 = 2ml (assuming equal quark masses). 



be more precise, one should distinguish random matrix parameters and physical parameters in the relations 2]Vzrmt = z = 
z phys 2 A/mRMT = m= m^y S VI., and 2N/i| MT = p. 2 = fi 2 ^ F 2 V. This distinction makes it explicit that the limits N — ■ oo and 
V — oo can be decoupled. 
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We now introduce the microscopic limits (denoted by a subscript s) of the orthogonal polynomials, 
the norm, and the weight function, respectively. They are worked out in App. A, and we obtain 



p v s [z; a) - lim 



Af-oo {2N) v+l12 
\ - 1™ f nAn2„2N„v 



p^_ 1 (z/2Af;d/2An = s/ne~ al2 z~ v I v (z) , 



r v s {&) = lim {2N)' L e iN r v N _ l {al2N) = 4n 2 ae & , 

AT— oo 

z 2 +z* 2 

w v Jz,z*;a) = lim (2N) 2v+2 w v [zl2N,z* I2N;&I2N) = |z| 2(v+1) e~ 8a K t 

N—oo I 4a 



(3.15) 
(3.16) 

(3.17) 



3.2.1 Quenched case 

In terms of the above definitions, the microscopic limit of the quenched average phase factor (3.7) is 
given by 



<e 2W ) = lim <e 2W ) 

\ K s /Nf=0 ~ flfc ^ /N f =0,al2N,ml2N 



rj(d) Jc z 



d 2 z 



■w*{z,z*;a) 



p v s (z*;a) z* 2 p v s +l {z*;a) 
p v s {fh;a) rh 2 p v s +l [rh; a) 



(3.18) 



As expected, the dependence on N has dropped out, leaving a finite microscopic limit for the average 
phase factor. Substituting the asymptotic results from Eqs. (3.15)-(3.17) yields 



(e 2W ) 



-2a 



IAz*) z*I v+l {z*) 
I v {rh) ml v+ \{m) 



4narh v Jc z 2 - rh 2 ' ' "\4a 
where we renamed the integration variable z back to z. This equation can be rewritten compactly as 

\ K s /N f =0 



(3.19) 



JCyfsia, rh) J€ v ,\(a,rh) 
I v {rh) ra/ v+ i(m) 



(3.20) 



where we defined the integral 



-2d 



Jfyjcia, rh) : 



4 



rylzl 'e 6a K v \ — \{z ) I v+ kiz ), 

4jiam v Jc z z - m z \4a 



(3.21) 



which is closely related to the microscopic limit of the Cauchy transform (2.12). For the quenched case, 
this integral is only needed for k = 0, 1, but as we shall see in the next subsection, in the unquenched case 
it will be needed for k — 0, . . ., Nf + 1. 

3.2.2 Unquenched case 

We now take the microscopic limit of Eqs. (3.12) and (3.14). In Eq. (3.12), the Vandermonde deter- 
minant Ajv /+ i(m 2 ,{m^}) is a product of Nf(Nf + l)/2 factors for which the microscopic limit yields a 
{2N) N f^ N f +1 ^ dependence on N, while the explicit mass and z* factors in the determinant yield a factor 
l/(2AD (Ar / +1)(iv / +2) . After also introducing the microscopic limits (3.15), (3.16), and (3.17) for the orthog- 
onal polynomials, their normalization factor, and the weight function, respectively, one finds 

N f l2 -alN f /2+2) 



Z, 



JV^+lll* 



{a,m;\m f }) = - [2 N) {2v+N f )N f 12 e~ NN f 



m 1 \ 4a 



4na{mrh\ rh2... mjv / ) v Ajv / +i {rh 2 , {rh 2 }) 

h.oiz*) I v ,\{z*) ■■■ I v ,N f +iiz*) 
I v ,o(fh) l v ,i(rh) ••• 7 VjW/+ i(m) 
Iv.oim) I v ,\{rhi) ••• / VjiV/+ i(mi) 



Iv,o(mN f ) Iv,i[rhN f ) ••• h,N f +\i.rh Nf ) 



(3.22) 
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where we have again renamed the integration variable from z back to z and introduced the notation 



(3.23) 



In the microscopic limit of Eq. (3.14), the Vandermonde determinant yields a factor (2N) { - N f~ 1)N f which 
exactly cancels the factor l/(2A0 (iV/ ~ 1)W/ coming from the explicit mass factors in the determinant. After 
introducing the microscopic limit (3.15) of the orthogonal polynomials we find 



n N f l2 e -&N f /2 



N, 



(mi m 2 ... m Nf ) v A Nf ({m 2 }) 



(3.24) 



with 



@v / ({%}) = 



Iv,o{rhi) I v ,i(rhi) 



Iv,N f -\{mi) 

Iv,Nf-li^2) 



Iv,o(rh Nf ) I Vl \{rh Nf ) ■■■ I ViNf -i{rh Nf ) 



(3.25) 



The microscopic limit of the average phase factor (3.9) is given by the ratio of Eqs. (3.22) and (3.24). The 
dependence on N drops out to give 



(e 2W ) 



-2a 



An&^YlJ^tfti - m 2 )Qiy f [{rhf]) 



(3.26) 



/ 



d l z 



-\z\ 



2(v+l) 



(z*y v K v 



4a 



Sa 



/v,0(z') Iv.liz*) 

I V fi{m) I v ,\{m) 



Iv,N f + l(Z*) 

Iv,N f +i(rti) 
Iv,N f +\{mi) 



Iv,o(m Nf ) I Vi i{m Nf ) ••• I ViNf+ i{m Nf ) 



which can be rewritten using the integral definition (3.21) as 

J€ v ${a,m) J€ v ±{a,m) 
Iv.oim) I v ,\{m) 
/v,o(mi) ly.iirhi) 



1 



U%irh 2 - m 2 )3>v f ({m f }) 



I v ,N f +l(m) 



(3.27) 



We have thus reduced the problem of calculating the phase factor to the calculation of the two-dimen- 
sional integral J€ y ^{a, m) in Eq. (3.21) for k — 0, . . .,Nf + 1. This integral will be computed in Sec. 4. 

3.2.3 Equal mass fermions 

We now consider Eq. (3.27) for the special case in which all dynamical fermions have the same mass m 
as the valence quark. To simplify the general expression we perform a Taylor expansion of the entries 
Iy^ifhf) of the determinant around m, 



Iv,kirhf) = I v ,k(m) + £ 



22 Mm 



(fiif- m) 



/= l,...,Nf. 



(3.28) 



Because a determinant remains unaltered when making linear combinations of its rows, we see that for 
each additional fermion it is sufficient to keep the next higher-order term in the expansion (3.28). The 
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lower-order terms will not contribute as they are identical to the contribution from one of the previous 
fermions in the determinant, while the higher-order terms can be neglected as their contribution will 
vanish when rhf — > m. After taking each fermion mass in turn to m, this leads to the simplified expression 



(e 2W ) 
\ e s IN, 



I v ,o(m) I v ,\{m) 

i' vo m r VA m 



{2m) N fN f \ 



i' V)0 m i' Vil m 

I 



Iv.Nf+iim) 



i lNf] m 



Iv,N f -l{m) 



I' 



(N f -l) 



m i v { m ••• i ' Am) 



(3.29) 



An alternative way to write this result is 



(e 2W ) = 

\ K s /N f 



Wjsr f {a, m) 



f {2m) Nf N f \W Nf [0,1,..., N f -l) 



(3.30) 



where we have defined 



Nf+l 



W N Aa,m)= £ {-) k je V}k {a,m)W Nf+l {Q,...,k-\,k+l,...,N f +l) 



k=0 



(3.31) 



as a sum of Wronskians of order Nf + 1 with indices ranging from to Nf + 1, where in each term a 
different index k is absent. The Wronskian 



W n [I Vikl {m),...,I Vikn {m)) ■ 



Iv^m) I v ,k 2 (m) ■■■ I v ,k n {m) 
I' . {m) I' . (m) ••• I' . (m) 

V,k\ v ' V,K2 v ' v,k n v ' 



I [n 7 l \m) I {n 7 l \m) ■■■ I [n 7 l \m) 

v,k\ v ' v,k-2 v,k n v ' 



(3.32) 



that appears in Eqs. (3.30) and (3.31) has been abbreviated by W n {k\,...,k„). 



4. Evaluation of the complex Cauchy transform 
4.1 Asymptotic behavior 

To investigate the two-dimensional integral (3.21) it is instructive to first study the asymptotic behavior 
of the integrand. For large values of its argument the TT-Bessel function behaves like [34, Eq. (9.7.2)] 




(4.1) 



Here and in the rest of the paper we use the ~ symbol for the leading-order term in an expansion for small 
or large argument, including all prefactors. To find the asymptotic behavior of the 7-Bessel function we 
first note that the modified Bessel functions satisfy the relation [35, Eq. (7.11.45)] 

K v {ze ±i7t ) = {-) v Ky{z) + inl v {z) . (4.2) 
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As the iT-Bessel function has a branch cut along the negative real axis, it is convenient to adopt the con- 
vention arg(z) e {—n,n\ for the complex variable z. According to this convention, reversing the sign of 
z = re' 6 yields 

[ re W-n) f or E (o,^], 

such that -z also has its argument in {-n,n]. Combining Eqs. (4.2) and (4.3) gives the relation 

I v {z) = l -^{{-) v K v {z)-K v {-z)) (4.4) 
n 

with 

{ + 1 for arg(z) £ (0,7r] , 
(4.5) 
-1 forarg(z)e(-7r,0]. 

Alternative definitions for rj{z) are ir){z) = \fz\\f^z or r\(z) - sgn(Imz), the latter only for z t R. Substi- 
tuting Eq. (4.1) in Eq. (4.4) gives the asymptotic formula 2 

iri(z) I .,e~ z e z \ 1 I e~ z e z \ 
ly(z) ~ (-) v — - — = — (-) v -= + -= , (4.6) 



271 \ \/~Z \T r Z) \fl~K\ V^-Z sfz 

where we used ir\{z) - \fzl\f zr z to derive the last expression. With Eqs. (4.1) and (4.6), the asymptotic 
behavior of the integrand in Eq. (3.21) is proportional to 



|z| 2v+1 




(4.7) 



where z = x + iy. Along the x-direction the integrand of the two-dimensional integral decreases like a 
Gaussian with width \f&. However, in the y-direction the integrand oscillates very rapidly inside an en- 
velope that goes like y v+fc - 3/2 . Eqs. (3.20) and (3.27) contain terms with v + k > 1, for which the integral 
(3.21) will diverge unless some particular cancellations occur due to the oscillatory behavior of the in- 
tegrand. As the integral represents an observable quantity in random matrix theory we do expect such 
cancellations to obtain a finite result. 

An instructive numerical exercise is the direct computation of the two-dimensional integral over 
the generalized Laguerre polynomials for finite N, as given in Eq. (3.7) for the quenched case. Although 
Mathematica can only handle the numerical integration for N < 30 because of the strong oscillations, the 
results for N - 1, . . . , 30 show a clear convergence towards a finite microscopic limit. (These numerical 
results also agree with the simulations presented in Sec. 5.5.) This clearly indicates that Eq. (3.19) is 
perfectly sane, even though the evaluation of the integral is nontrivial. 

4.2 Change of integration path and transformation of variables 

The main problem is to find a way to integrate over the oscillatory behavior in the y-direction. For v = 
a method was devised in Ref. [28] in which the integration along the real y-axis in the original (x, y)- 
plane was deformed to an integration path in the complex y-plane. This results in a well-behaved one- 
dimensional integral. Here, we show how this derivation can be generalized to v ^ 0, where proper care 
has to be taken of an additional singularity occurring in the integration domain. 

The two terms e~ z and e z in Eq. (4.7) behave differently for | y| —> oo. The first term decreases expo- 
nentially in the upper half of the complex y-plane and diverges exponentially in the lower half, whereas 
the second term behaves the other way around. Rather than treating these two terms separately, we can 



2 Note that the asymptotic formula (9.7.1) in Ref. [34] only contains the second term of Eq. (4.6) and cannot be used here as it is 
only valid for | argz| < it 12. 
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+£ 



y r 



Figure 1: Deformation of the y-integration from the original integration along the y r -axis to a path on which the 
integrand vanishes sufficiently rapidly for |y| — > oo. 



use Eq. (4.4) to simplify Eq. (3.21). Because of the z — ► -z symmetry of the integrand in (3.21), the two 
terms in (4.4) give the same contribution to the integral, and we obtain 



JCvjcia, m) = ■ 



ie~ 2& f 
2n 2 arh v Jc 



\z\' 



We now deform the y-integration path from the real y-axis to the path shown in Fig. 1, where y = y r + iyt 
and £ — ► after the integration over d 2 z. This can be done since the integrand vanishes exponentially 
for | y\ — * co along the deformed path and since the integrand has no singularities between the real y-axis 
and the deformed path. Writing y = is + e on parts A and C of the path, respectively, we find for Eq. (4.8) 

/oo I rO nE poo \ 

dx\i \ ds f{x, is- e)+ / dy r f(x,y r ) + i \ dsf(x,is + e)> 
■oo I Joo J-£ JO ) 

/oo poo 
dx I ds [f+(x, is) -/_(jc, is)] , 
■oo JO 



(4.9) 



where we introduced the notation f+ {x, is) — f{x,is± e) with 



fix, is) - 



ie -2& ( _ ) v+k (x _ s] v+i 
2n 2 am v ix-s) 2 -m 2 



4& K v 



4a 



r]ix+ s) ix+ s) k+1 K v+k ix+ s) 



(4.10) 



Note that when continuing y to the complex plane we have rewritten the integrand as an explicit function 
of x and y, since z and z* are no longer complex conjugate. The second integral in the first line of (4.9) 
gives zero in the e —> limit since the integrand is regular at y = 0, which is most easily seen by looking 
back at Eq. (3.21). As e — «■ in Eq. (4.9) we are thus left with the difference of two integrals over semi- 
infinite sheets infinitesimally close and parallel to the (x, s) -plane, which we denote by S+ and S_ for 
y r > and y r < 0, respectively. 

The integrand (4.10) can be further simplified by introducing the variable transformation 



t - x- s 
u- x+ s 



x=(u+ f)/2 
s-iu- t)/2 



(4.11) 



with Jacobian 1/2, yielding 

/oo roo 
dt I du[f+it,u) + f-it,u)\, (4.12) 
■oo it 

where the integration limits in the transformed variables can be read off from Fig. 2. We introduced the 
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Figure 2: Structure of the integration plane. The integral is computed over the semi- infinite sheets S±, with s > 0, 
just above and below the (x, s) -plane. The shaded dot labeled by A at m = 0, t = -m indicates the singularity of the 
integrand responsible for a new contribution to the phase factor that arises when v ^ 0. 

notation f±[t, u) = f[t+, u+) with t+ = x- s+ ie and u+ — x+ s+ ie, corresponding to the integrand on the 
two sheets S+. In the transformed variables the integrand is given by 

e -2a ,y+k t v+l u k+l g+jgj ( fu\ 

f{t,u)= — — e 8& Kv — \K v+k (u), (4.13) 

where we have absorbed the Jacobian and factors of i in the definition of /. We have used T]{u+) = +1 
according to Eq. (4.5), which changes the difference in Eq. (4.9) to a sum in Eq. (4.12). 

4.3 Singularities of the integrand 

When taking the limit £ — ► in Eq. (4.12) one has to take into account the singularities of the integrand in 
the [t, z/)-plane. There are two mass-pole lines parallel to the w-axis at t — +m as well as singularities and 
branch cuts of the iC-Bessel functions for zero and negative real argument, respectively (see Fig. 2). 

For future use we first analyze the branch cut of the Bessel functions along the negative real axis. 
When going from one Riemann sheet to the next across the branch cut, K v {z) changes by {-) v+l 2nil v {z) 
[35, Eq. (7.11.45)]. For ieI we have 

\K v (\x\) forx>0, 
lim K v {x±ie) = \ (4.14) 
£-o+ \{-) v K v {\x\) + inI v {\x\) forx<0. 

Consider separately the three sectors of the integration region indicated in Fig. 2. Using (4.14) together 
with t± - t+ ie, u± = u+ ie, and {tu)+ — tu± ie (with e > because s > 0), we can write the product 
K v (tu)K v+ k(u) (with the factor l/4a omitted for simplicity) appearing in (4.13) as 



lim [K v {tu)K v+ ic(u)] + = < 

£—0+ 



K v {\tu\)K v+k {\u\) foit,u>0, 
(-) v K v (\tu\)K v+k (\u\) + i7iI v {\tu\)K v+k (\u\) fort<0<u, (4.15) 
(-) v+k K v {\tu\)K v+k {\u\) + inK v {\tu\)I v+k {\u\) for t,u<0. 



Note that the fourth quadrant, where t > and u < 0, does not overlap with the region of integration and 
thus does not need to be considered. 
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-A K±J s+ 

C + 

Figure 3: Deformed integration paths over u for t= —m. 



We now turn to the mass-pole factor in Eq. (4.13), which can be written as 

1 1 ( 1 1 

{t+ie) 2 -rh 2 2m\t+ie-m t+ie+m 

This enables us to apply the Sokhatsky- Weierstrass theorem 

rb fi x \ rb fr x i 

lim dx=+inf(0) + PV I - — dx, (4.17) 

E-0+J a X+ie J a X 

where PV denotes the Cauchy principal value integral, to perform the r-integral in Eq. (4.12), which thus 
becomes a sum of the residues at the mass poles t-+m and a principal value integral over the complete 
f-axis. In App. B we show that the principal value part of the t -integral in Eq. (4.12) vanishes because of 
symmetry properties so that the r-integral is entirely determined by the mass-pole contributions. The 
residue contributions from the mass poles yield two one-dimensional integrals over u along the f = + m 
lines, resulting in 

in \ f°° f°° } 

.X v k {a,m)- { lim I du\g(jfi, u+) —g(fh, U-)\ — lim / du\g(—m,u + ) — g[—m,u-)]>, (4.18) 

2m U-0+ Jm e^O+J-m J 

where 

git, U) = {t 2 - m 2 )f{t, u) = --^—[-) v+k t v+1 u k+1 e- L ^~ K v f^l K v+k {u) . (4.19) 

While the first integral in Eq. (4.18) is well-behaved, for v > the integrand of the second integral (for 
which t = - m) is singular at u = 0. This singularity is labeled by A in Fig. 2. To avoid this singularity when 
taking the £ —> limit, we deform the integration path on S± near zero as shown in Fig. 3 and then take 
£ — ► 0. This yields 

in \ r°° r r~ s r°°i ] 

^f VJ t(d,m) = < | duG{m,u) - lim J + | \duG{-m, u) \ + A v ^{a, m) (4.20) 

2m [Jm 5-0+LJ-m Js J ) 

with G(f, u) = lim £ ^o + lg(t, u+) - g{t, U-)}. The contribution of the singularity at u — 0, which for simplic- 
ity will be called the A-term, is 

A VJ t(d,m) = lim \ I dug{-m,u)- \ dug{-m,u)\ , (4.21) 

2 m 5-o+ ( Jc + Jc_ ) 

where C± denotes the two small semicircles shown in Fig. 3. We shall see in Sec. 4.5 that this term makes 
a contribution to the phase factor for v ^ 0. 

4.4 Contribution of the branch cut discontinuity 

The integrand in the curly braces of Eq. (4.20) involves differences that can be simplified using Eq. (4.15), 

for t, u > , 

-2inI v (\tu\)K v+k {\u\) forr<0<ii, (4.22) 
+2inK v {\tu\)I v+k {\u\) for t,u<0, 



(4.16) 
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where we have again omitted the factor 114a. The only contributions in Eq. (4.22) come from the branch 
cut discontinuity of the iC-Bessel function for negative real arguments. The first case in Eq. (4.22) vanishes 
as neither Bessel function has a negative argument and no branch cut discontinuity is encountered. We 
now use Fig. 2 to distinguish the various integration regions of Eq. (4.20) . For the t-m integral we always 
have t, u > 0, and according to Eq. (4.22) this integral vanishes. Substituting Eq. (4.22) in Eq. (4.20) with 
integrand (4.19), we are left with 

. 2 

„ e ~ 2a -M [ f» k , , (\mu\ 

je Vik {a,m)= — 7T — |J rfu(-) u e ««/, 



4a [Jo [ 4a 



+ J du{-) k+1 u k+1 e 8&K v ^^jI v+k (\u\)^ + A v ,ki&,m, (4.23) 

where the integration limits +5 in the curly braces of Eq. (4.20) were set to zero as the integrands are 
regular for u — ► 0. This expression can be simplified by transforming u — ► - u in the second integral, and 
we obtain 

. 2 

e~ 2 "~w I r°° i , , _m*l I mu\ C" 1 t , _!£: I mu\ ) 

je Vik {a,m) = — — — |J du{-) k u k+v e a& j v I — I K v+k {u) + J duu k+v e a&K v I — J I v+k {u)\ 

+ A v , k {a,rh). (4.24) 

For v = this expression is identical to the result given previously in Ref. [28], as the A-term vanishes in 
this case. For v ^ the additional contribution is important and will be computed below, see Eq. (4.33) 
for the result. We shall see that it even dominates for small m. 

4.5 Contribution of the Bessel function singularity 

The semicircles around the singularity in Eq. (4.21) run in opposite directions on the sheets S+ and S_. 
Reversing the direction of integration on one of the sheets changes the difference in Eq. (4.21) to a sum. 

The JT-Bessel functions in the integrand can be split into a meromorphic part with a pole at u — and 
a part containing the branch cut. For 5^0 the integrals along C+ would not only receive contributions 
from the singularity at u = but also from the branch cuts, lying on both sides of the origin, which were 
already included in the line integrals given in Eq. (4.24). When 8^0 the branch cut contributions will 
vanish from the integral (4.21), and only the meromorphic part of the integrand will contribute to the 
A-term. Therefore we can rewrite the A-term (4.21) as an integral over a closed contour Fo enclosing the 
singularity of the meromorphic part g(-m, u) of the integrand. We thus obtain 

in f 

k vk {a,fh) = <b dug{-m,u), (4.25) 

2m h a 

where Fo consists of C+ and C_ and is traversed in the counterclockwise direction. 

For z — ► 0, K v {z) diverges logarithmically for v = and like z~ v for v > 0. Upon multiplication by 
u k+l , for v = the integrand (4.19) is not singular at u — 0, and thus the A-term vanishes in this case. For 
v > the integrand has a pole of order (2v - 1) at u — 0. Using the residue theorem, we obtain 

in n 2 

A vk (a,m)= 27r/ Res g= a\ , (4.26) 

2m u=o m 

where a\ is the coefficient of the u~ 1 term in the Laurent expansion of g(- m, u) around zero. To find this 
coefficient, we neglect the ^-independent terms in Eq. (4.19) and write the ^-dependent terms as 

/ u 2 \ t-mu\\ p 1 
uf(u)g(u)Hu) = uexp[-— ) k v[-^t-) [u lc K v+k (u)\ . (4.27) 
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As the product ug[u)h[u) has a singularity of order 2v - 1, we need to perform the Taylor expansion of 
the exponential to order 2v - 2 to find all contributions to the simple pole. Hence 



8a j £? £l { 8a 
For v > the series expansion of K v {z) for small z is [34, Eq. (9.6.11)] 

2 v - lv - 1 (v-l-fc)!f z 2 \ k ^ v 
W = — E— (- T ) + (4-29) 

The functions g and /i containing the Bessel functions each need to be expanded to order v - 2 using 
Eq. (4.29), as the remaining part of the integrand is as singular as u 1 ~ v . So 

2 V " 1 v "! (v-l-i)! f (&w) 2 V 

g(u) s K v {bu) = — — £ — + 0{u v ) (4.30) 

[bu) v f^ Q i! I 4 j 

with = - rhIAa and 



h{u)=u k K v+k {u)= — V n —7" +<?(")■ (4-3D 

" v jTo J- \ A 



Putting things together yields a triple sum, given by the product of Eqs. (4.28), (4.30), and (4.31). Since we 
are only looking for the coefficient of u~ 1 , we obtain the condition 1 + 2( + (2 i - v) + (2j - v) = - 1, which 
eliminates one of the sums by setting £ = v - 1 - i - ;'. Since £ > we also have i + j < v - 1 and thus 

Resl ufiu)g(u)Hu)} = - !+ y V 1 tv - 1- 0!(v + fc - 1 -7)! 2V+fc+1 _ 3; . +j _ v+2i - _ ; . +J - +1 ^ 
"=° ij^o (v-l-i-;)!/!;! 

Combining this expression with the (/-independent terms in Eq. (4.19) and substituting it in Eq. (4.26) 
yields 

M * 2 (-) fc 2 v+fc " 1 ! ' + ^ v " 1 ( v -l-i)!(v+jfc-l-/)! /m 2 V 

A v , fc (a,m) = e - 2a -w^_^ £ . . J Upr (2 a)', (4.33) 

m v (v-l-i -;)!;!;! UaJ 

which is a bivariate polynomial of degree v - 1 in m 2 la and a which contains v(v + 1) 12 terms. Our final 
result for J€ v ^{a, rii) is thus 



-2a " 



■^f vfc (a,m) = if du(-) « e 8i/J U: v+fc (u) + / duu k+i e a&K v \ \l v+ 

4d I Jo I 4d J Jo l4aj 



fc(") 



+ e 



2 „_f (-j^^- 1 '^f - 1 (v- 1 - iW+k-1 - ;)! I m 2 , , 



i,;=o 



(v- 1- i-j)W.jl [8& 



Y — — — {2&) ] . (4.34) 



As mentioned earlier, for v < one just needs to replace v by |v|. Our general expression can also be used 
for v = as the double sum vanishes and the correct v = result is reproduced. 

Note that for k = 0, Eq. (4.34) has been computed in Ref. [36] using a different method. The re- 
sult given in Eq. (65) of that reference looks rather different from our result but agrees numerically with 
Eq. (4.34) for k — after adjusting some prefactors. However, for the calculation of the phase factor using 
Eqs. (3.20) and (3.27) we need the more general result of Eq. (4.34) for arbitrary k. The method of Ref. [36] 
does not straightforwardly extend to k ^ since it uses orthogonality relations that only hold for k — 0. 
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5. Explicit results 
5.1 Quenched case 

We now write down the explicit form of the phase factor for the quenched case using the solution for the 
Cauchy transform integral derived in the previous section. Recalling Eq. (3.20), we obtain 



(e liB ) 

\ e s /Nf=0 



-2a- 



4a 



I du 

Jo 



ue 8d l 



mu 
4a 



~8&K V 



mu 
4a 



ml v+ i {m)K v {u) + I v (m) uK v+ i{u)] 
[ml v+ i {m) I v [u) - l v (m) u/ v +i(u)] \ + A 



where 



N f =0, 



A, 



[a, in) — e 



J duue 8( 
Jo 

2& _nl 2 v ~ l (v- 1 - i)!(v- 1 - ;)! I m 2 ) 1 



{a, m) , 



E 

i.j=0 



{v-l-i-j)liljl 



8a 



1 ~2{v-j) 
I v [m) mly + i [m) 



(5.1) 



(5.2) 



The contribution of both terms to the total phase factor will be illustrated in the numerical results of 
Sec. 5.5. 

5.2 One and two dynamical flavors with equal masses 

We now give explicit expressions for the one- and two-flavor case. For one dynamical fermion with mass 
equal to the valence quark mass, the phase factor from Eq. (3.29) is given by 

J£yfi[a, m] J€ Vi \{a,m) J€ Vt 2{a,rh) 
7 v>0 (m) I v ,i{m) I v , 2 {m) . (5.3) 

K.o^ K.iW K,2^ 
Substituting the solution (4.34) derived in Sec. 4 for the complex integral we find 

K v {u) -uK v+ \{u) u 2 K v+ 2{u) 
I v ,o(m) l v ,\{m) I v ,2{m) (5.4) 

i' v0 m i' Vtl m i' v>2 m 



(e 2W ) 

\ K 5 /Nf=l 



1 

2ml v {m) 



-2a- 



(e 2W ) 

{ s }N f =1 8am 



l ~m I r c 
Ivirh) | Jo 



duue 8d/ 



mu 
4a 



SaK v 



mu 
4a 



+ I duue Be 
Jo 



Iy(U) Uly + \{U) U 2 Iy +2 {U) 

Iy.oim) ly,i{m) I Vi 2(m) 

K.o^ K.i^ K,2^ 

1 -2(v-j) 



m v+1 Iy{m) 



E 

i,j=0 



(v-1- i- jl 



8a 



(2a) J 



Iyflim) I v ,i(m) 

I'm I'm 



4{v-j) 2 
I v ,2{rh) 

I'm) 



with the Pochhammer symbol (a) n defined in Eq. (C.18). Using I' v {z) - I v -\{z)-vl v {z)l z [34, Eq. (9.6.26)] 
the derivatives in Eq. (5.4) can be explicitly computed, 

I' vk m= [m k Iy +k {m)}' 



m k ~ l y m/ v+fc _i (m) - vl v+k {m) 
For two dynamical fermions with masses equal to that of the valence quark, Eq. (3.29) yields 

J?vfi[a, rh) Jtfy t \[a,rh) J€ Vt 2{a,rh) J€ v ${a,m) 



(5.5) 



(e 2W ) 

\ K s /N f = 



I v ,o(m) 

i' Vi0 m 
i" m 



I v ,i(m) 

i' Vtl m 
i" A m 



I v ,2{m) 

i' Vt2 m 
rum 



Iy,3m 

i'y, 3 m 

I'Um) 



Iy.om Iy.lim) 

i' v0 m i' Vyl m 



(5.6) 
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One can show that 



7 v ,o(m) 7 v ,i(m) 

= m\I v {m) -7 v _i(m)7 v+ i(m) , 

7; (m) 7; ;1 (m) 

and using the solution (4.34) for the complex integral we find 



(5.7) 
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m v+s \l v {m) 2 - 



'^iT 1 (v 

i(m)7 v+ i(m)] ( -j =0 



(v-1- 

1 -2(v-j) 4(v-y) 2 
7 v , (m) 7 v ,i(m) 7 Vj2 (m) 
7; (m) 7; x (m) l' v2 {m) 
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rum 



j)\ I m 2 
8a 



(2a) 



(5.8) 



with the Pochhammer symbol [a) n from Eq. (C.18). 

In Fig. 4 we compare the mass dependence of the average phase factor in the quenched case, given 
byEq. (5.1), with the predictions for one and two dynamical flavors from Eqs. (5.4) and (5.8). Although the 
sign problem becomes less severe as the fermion mass increases, the dynamical quarks have a negative 
effect as they clearly reduce the value of the phase factor. Note that, surprisingly, for v = and small mass 
[m < 1) the effect of the dynamical quarks seems to be reversed (this is true for any value of a). This is 
related to the fact that for v = 1 the chiral limit of the average phase factor is equal to e~ 2a , independently 
of the number of flavors, see Eq. (5.10) below. This is also verified in the middle plot of Fig. 4. 



< c f>> 
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N f =0 


v=0 




















Figure 4: Average phase factor of the fermion determinant for Nf = 0, 1,2 as a function of the fermion mass m for 
d = 1.0 and v = 0,1,2. We have verified that the m -» limit of the curves agrees with the chiral limit computed 
analytically in Sec. 5.3. 
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5.3 Chiral limit 

In the limit m — ► the phase factor (3.30) can be simplified and written in terms of special functions. The 
limit has to be derived differently for v = and for v ^ 0, and the general derivation for arbitrary Nf can 
be found in App. C. For trivial topology the result is given by Eq. (C.15), 

< e f >v=o,m=o = Wf + Di2a) N f +1 T{-N f - 1,2a) , (5.9) 
where T{a, z) is the incomplete gamma function. For nontrivial topology the result is given by Eq. (C.17), 

< e ? >v>0,m=0 = e E —f^ " 1} * (5.10) 

where the Pochhammer symbol {a) n is defined in Eq. (C.18). Note that for v = 1 the phase factor is 
simply given by e~ 2a , independently of Nf. The result (5.10) can also be expressed in terms of incomplete 
gamma functions as shown in Eq. (C.25), 



Nf + l 



1 'V 

{ s )m = (v + N f V. t. 



k 



(v-k) Nf+l (2a) k Y(v-k,2a). (5.11) 



The last expression is also valid for v = 0, see Eq. (C.26). The quenched results are obtained by setting 
Nf - in the above equations. 

As discussed in App. C, for m = and v ^ the phase of the determinant is exclusively given by the 
new A-term. It is interesting to note that the contributions to the phase of the determinant in the chiral 
limit originate from different terms for v = and v ^ 0. 

5.4 Thermodynamic limit 

The thermodynamic limit of the phase factor (3.30) is the limit of that equation for a = fL 2 F 2 V — » oo and 
m — mVL — ► oo. As discussed in Ref. [28], this limit depends on whether 2a/m (or, equivalently 2^lm n ) 
is smaller or larger than 1. We find that for 2a < m the thermodynamic limit is given by 

< e f e ) th =(l-2a/m) N / +1 . (5.12) 

As expected, the thermodynamic limit does not depend on v, which is a consistency check of our result. 
Equation (5.12) agrees with the special cases v = and Nf = 0, 1,2 considered in Ref. [28]. The proof of our 
general result is given in App. D. Note that the contribution of the A-term vanishes in the thermodynamic 
limit. For 2a > m the phase factor is exponentially suppressed in the volume so that its thermodynamic 
limit is zero. An asymptotic large-volume expansion of the phase factor could be computed for this case 
from Eq. (3.30), analogously to Ref. [28]. 

In Fig. 5 we show how the thermodynamic limit is approached for the case of Nf = 2 and v -2. 




Figure 5: The left plot shows how the a-dependence of the phase factor changes as a function of m. In the right 
plot the roles of a and m are reversed. The curves m — > oo and a — oo correspond to the thermodynamic limit 
(5.12). The dashed lines indicate the contribution of the A-term to the phase factor (this contribution vanishes in 
the thermodynamic limit). 
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5.5 Numerical simulations 



To check our analytical results and especially the contribution of the A-term, which is new in the v ^ 
case, we performed numerical random matrix simulations to compute the average phase factor of the 
fermion determinant in quenched chRMT. The simulation details can be found in App. E. 

In order to keep the statistical error sufficiently small the simulations were performed with samples 
of 100,000 random matrices. Figure 6 shows the a-dependence of the quenched average phase factor 
for v = 0, 1,2,3,4,5 with m = and matrix size N = 20. The simulation results are compared with the 
chiral limit of the analytical results given in Eqs. (5.9) and (5.10) with Nf = 0. The agreement is extremely 
good except for larger values of v (> 3) where the simulation results lie slightly, but systematically, below 
the predictions. This is merely a finite-Af effect, which is also clear from the comparison with the exact 
finite-A? result shown for v = 5. 




Figure 6: Average phase factor of the fermion determinant with to = in the quenched case for varying chemical 
potential parameter a and v = 0, 1,2,3,4,5. The full lines are the predictions of Eqs. (5.9) and (5.10) for Nf = 0. For 
v = 5 we also show the exact finite-iV result from Eq. (3.7). The data points were computed from RMT simulations 
with matrix size N = 20 and 100,000 samples. No error bars are shown since they are smaller than the data points. 



The convergence towards the microscopic limit is il- 
lustrated in Fig. 7, where we show the Af-dependence of 
the average phase factor for two typical cases: fast con- 
vergence for a - 0.22, v = 2 versus slow convergence 
for a = 2, v = 4. The figure can help us determine how 
large the random matrices need to be in order to repro- 
duce the analytical results in the N — «■ oo limit. In the fig- 
ure we also show the A/-dependence of the phase factor 
from the theoretical framework, by numerically solving 
the two-dimensional integral for the finite -N expression 
(3.7), which is expressed in terms of generalized Taguerre 
polynomials. We find perfect agreement with the data of 
Fig. 7, within statistical errors, for N from 1 to about 30, at 
which point the integrals oscillate too strongly, prohibit- 
ing Mathematica from performing the numerical integra- 
tion with sufficient accuracy. 



<<£"> 

0.8 



a=0.22, v=2 








a=2.0, v=4 





10 20 30 40 JV 50 

Figure 7: Average phase factor vs matrix size N. 
The horizontal lines show the analytical predic- 
tions in the microscopic limit (AT — > oo), while the 
solid lines going through the data points are the 
finite-iV results from Eq. (3.7) for N < 30. The up- 
per curve was computed with a = 0.22, v = 2, the 
lower with a = 2, v = 4, both for Nf = and m = 0. 
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Figure 8: Average phase factor of the fermion determinant as a function of the chemical potential parameter a for 
v = 3, 4, 5 with Nf = and m = as in Fig. 6, but with increased N = 80. 



From the study of the N- dependence we conclude that for v = 0, 1, 2 it suffices to take N - 20. How- 
ever, for v = 3, 4, 5 we increased the matrix size to N = 80 to be close enough to the microscopic limit. As 
can be seen in Fig. 8, this increase in N results in excellent agreement between simulations and analytical 
predictions for larger v as well. Note that the m = case is a key test for the new v ^ contribution (5.2), 
as only this term contributes in the chiral limit. 

From Figs. 6 and 8 we conclude that, as expected, the average phase factor becomes unity when the 
chemical potential vanishes, as the Dirac operator then becomes anti-Hermitian and the determinant 
real. For large a the average phase factor goes to zero, pointing to the increasing sign problem in dynam- 
ical simulations at large chemical potential. Observe that for increasing topology the sign problem seems 
to be delayed, as it sets in at a larger value of the chemical potential. 

In Fig. 9 we verify the mass-dependence of the analytical formula (5.1) and compare its predictions 
with the results from random matrix simulations as a function of m for fixed a. Note that the convergence 
to the microscopic limit slows down as the mass increases. This is noticeable in Fig. 9, where the N — 20 
data (red points) show a systematic deviation from the RMT-predictions. When increasing the size to 
N = 80 (blue squares) the agreement improves substantially. The importance of the new A-term for v ^ 




2 4 6 8 f n 10 2 4 6 8 f n 10 2 4 6 8 10 



Figure 9: Average phase factor of the fermion determinant as a function of the fermion mass m for Nf = 0, a = 1.0, 
and v = 0, 1, 2, 3, 4, 5. The simulations were performed with N = 20 (red points) and N = 80 (blue squares) . The latter 
data are already very close to the RMT-predictions for N — > oo (full lines). The filled area shows the contribution 
of the A-term specific to v 5^ 0. (Note that we generated different random samples for each m to avoid misleading 
correlations between measurements.) 
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is highlighted in Fig. 9 by the gray area, which corresponds to the contribution of the A-term in Eq. (5.1). 
This clearly shows how this term dominates for small masses. We also observe in Fig. 9 that for fixed m 
and & the sign problem becomes less severe as the topological charge is increased. 

Figures 6-9 demonstrate that the numerical simulations confirm the analytical prediction (5.1) for 
general topology in the quenched case. We are currently investigating the implementation of unquen- 
ched random matrix simulations at nonzero chemical potential using the analytical information we have 
obtained about the sign problem. 



6. Conclusions 

Dynamical lattice simulations of QCD at nonzero baryon density are plagued by the sign problem caused 
by the oscillating fermion determinant. To investigate this problem it is helpful to employ the equiva- 
lence between chiral perturbation theory in the e-regime of QCD and chiral random matrix theory, which 
also holds at nonzero chemical potential. As the average phase factor of the fermion determinant is an 
important clue in the study of the sign problem, it is a valuable quantity to compute in the framework of 
chiral random matrix theory. 

In this paper we derived an analytical formula for the average phase factor of the fermion determi- 
nant in quenched and unquenched chiral random matrix theory for general topology. The formula is a 
nontrivial extension of the result previously published by Splittorff and Verbaarschot for zero topology 
[28]. For nonzero topology a new contribution shows up, which dominates the phase factor for small 
valence quark mass. The new formula suggests that the severity of the sign problem is reduced as the 
topological charge increases. We also computed the chiral and thermodynamic limits from our general 
formula. 

The quenched formula was verified by random matrix simulations in different regimes of mass and 
chemical potential and for different values of the topological charge. Excellent agreement was found be- 
tween theory and simulations. We are currently in the process of comparing the RMT predictions derived 
in this paper to lattice QCD data computed with the overlap operator at nonzero chemical potential. 
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A. Microscopic limit of the orthogonal polynomials 

In this section we define and compute the microscopic limits of the orthogonal polynomials (2.9), the 
normalization factor (2.11), and the weight function (2.7). 

The microscopic limit of the orthogonal polynomials (2.9) is defined as 

e N e N I a l^" 1 (JV-1) 1 ( z 2 \ 

p v Az;a) = lim —p v N Azl2N;al2N) = lim — 1 = V N A . 

Fs Af-oo (2iV) v+1/2 iv-oo (2AT) V+1/2 I 2NI N N ~i w_1 \ 4N j 

(A.1) 

Using the definition of the exponential function, 

lim (l-— ) =e~ &12 , (A.2) 
W-ool 2N) 



-22- 



and Stirling's formula, from which we obtain 



we find 



™\f2nNe~ N , (A.3) 



z 2 



p v s {z;a) = lim ^e~ a l2 (2Ny v L v N _ 1 -— . (A.4) 

oo \ 4jV / 

The N — «■ oo limit of the Laguerre polynomial in Eq. (A.4) is given by [37, Eq. (8.982.2)] 

lim N- v L v J-^-\ = 2 v z- v / v (z) , (A.5) 
iV— *oo JV ^ 4 AT; 

and hence Eq. (A.4) becomes 

p v s (z; a) = ^fne~ &l2 z~ v I v (z) . (A.6) 
Next, the microscopic limit of the normalization factor (2.11) is defined as 

(2N) 2 e 2N a I a \ 2N ~ 2+V 
r v Aa)= lim {2N) 2 e 2N r v N , [&I2N) — lim - — i- n — 1 + — (N- \)\{N - 1 + v)! 

S AT-oo N ~ l AT-oo N 2N+V 2 N { 2N j 

= 4n 2 ae & , (A.7) 

where we again used Eqs. (A.2) and (A.3). Finally, the microscopic limit of the weight function (2.7) is 
defined as 



2v+2,,.v, s/ 

2 , **2\ / Am i a, /I An i*|2 



w v Jz,z*;a) = lim [2N] w v {zl2N,z* l2N;al2N) 

AT— oo 



lim l^i'cxpf Ml-fl/2JV)^ + r^ f N(l + a/2JV) |z|M 
tf-oo F i 4a/2AT 4N 2 J v i 2d/2AT 4N 2 ) 



N- 

'^ +1) exp|-^j^|-|. (A.8) 



z 2 + z* 2 \ f|z| 2 



B. Principal value integral 

In this section we show that the principal value integral originating from the application of the Sokhatsky- 
Weierstrass theorem (4.17) to Eq. (4.12) vanishes because of symmetry considerations. 

For v > 0, the integrand (4.13) of the principal value integral over t is singular along the line u — (for 
which t < in the integration region). Therefore we split the M-integral into a principal value part and a 
line integral circumventing the singularity of the Bessel functions at u = 0. (To simplify the notation we 
do this also for v = even though there is no singularity in this case.) For any t < 0, the w-integration in 
Eq. (4.12) is thus rewritten as 

noo I n—d f r°°] 

J duf±[t, u) = lirn jj +J + J ^duf±(t,u), (B.l) 

where C± denotes the small semicircles shown in Fig. 3. The principal value integral over t in Eq. (4.12) 
therefore becomes 

/oo roo roo roo 

dtliml du\f + {t,u) + /_(£, m)] = PV tu / dt I duF{t,u) (B.2) 
-oo 0+Jt J-oo Jt 

+ PV f f dt lim | f duf+{t,u)+ f duf-{t,u)\, 
J-oo 5-0+ I J c+ Jc_ ) 
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where F[t, u) = lim £ „o+ [/+(£> u) + f-(t, u)\. The total principal value integral W tu can be split, for any t, 
as 

/oo poo poo poo poo pO 

dtj duF{t,u) = VV tu \ dtj duF{t,u) + PV tu I dt I duF{t,u) . (B.3) 
-oo Jt J—oo JO J—oo Jt 

" v ' * v ' 

= A =B 

From Eq. (4.13) we see that the dependence of the integrand on t and u is given by 

t v+l u k+l t 2 + u 2 I tU \ 

f[t,u)<x— ? — 8d K v — \K v+k {u), (B.4) 

and according to Eq. (4.15) the sum of the contributions on the upper and lower sheets, S+ and S_, is 
proportional to (again with the factor l/4a omitted) 



lim [K v (tu)K v+k (.u)].+ [K v (tu)K v+k {u)]_ 



2K v {\tu\)K v+k {\u\) for t,u>0, 

2{-) v K v {\tu\)K v+k {\u\) forf<0<M, (B.5) 
2{-) v+k K v (\ tu\)K v+k (\ u\) for t,u<0. 



From Eqs. (B.4) and (B.5) we see that the integrand of A in Eq. (B.3) is odd in t, and hence the principal 
value of that integral vanishes. Integral B of Eq. (B.3) can be rewritten as 

dt I duF{t,u)-VW tu \\ dt I duF{t,u)- I dt I duF{t,u)\ 

■oo Jt l J-oo Jt Jo Jo ' 

Uoo rt poo rt -. 

dtj duF(-t,-u)-J dtj duF{t,u)Y (B.6) 

From Eqs. (B.4) and (B.5) we find for u, t > 

F{-t,-u) = {-) v+ \-) k+1 {-) v+k F{t,u) = F{t,u), (B.7) 

and hence (B.6) and the total principal value integral (B.3) vanish because of symmetry considerations. 

Next we treat the last term in Eq. (B.2). For v = the integrand (B.4) is regular at u — 0, and thus the 
integral trivially vanishes as 5 — «■ 0. For v > one can expand the integrand in pole terms behaving like 

u~ 2{v ~ i)+l , / = 0,...,v-l (B.8) 

for u — «■ by using the small- argument expansion (4.29) of the TT-Bessel functions. The w-integration 
around the pole flows in opposite directions for the integration on the lower and upper sheet. The simple 
pole gives contributions with opposite signs on the two sheets so that their sum is zero. The higher-order 
poles of the integrand are also odd in u, and one can show that the integrals along C+ and C_ vanish 
individually, so that 

PVj f dttimlf duf+{t,u)+ f duf-{t,u)\ = Q. (B.9) 

J-oo 5-0+ ( Jc + JC- J 

This completes the proof that the integral (B.2) vanishes. 



C. Chiral limit 

In the limit m — ► the phase factor can be simplified and written in terms of special functions. Our 
starting point is Eq. (3.30), and we need to compute the chiral limit of the Wronskian (3.32), which con- 
tains the function I Vik defined in Eq. (3.23). For small argument m the leading order term of the 7-Bessel 
function is [34, Eq. (9.6.10)] 

m 

I v {m) , (C.l) 

v ' 2 v v! 
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from which we obtain 



and its p-th derivative 



m v+2k 



I vk m ; (C.2) 

2 v+k {v+k)\ 



M m - - (v + 2fc)! fn^-v . (c.3) 

v,k 2 v+fc (v+fc)!(v + 2fc-p)! 

Substituting these expressions in the Wronskian (3.32) and using properties of the determinant gives the 
leading- order result 

y^nv+2Y.\ ki-n{n-\)l2 

where A n (ki,...,k n ) is a Vandermonde determinant. From this we compute the chiral limit of the de- 
nominator in Eq. (3.30), 



W N A0,l,...,N f -l) rr-j n £ -' ^ C - 5 ) 



2 JV / v n;: / /=1 

where we used the identities 

IVf-l JV/-1 

A JV/ (0,l,...,A/' / -l)= f] ^! and £ i = -N f (N f - 1) . (C.6) 

^=1 i=0 2 

From Eq. (C.4) it is easy to see that in the limit m — ► only the Wronskian corresponding to the k-Nf+1 
term in Eq. (3.31) will contribute to leading order, while all other terms will be of higher order. The phase 
factor (3.30) can therefore be written as 



( _)JV /+ i W Nf+1 {0,l,...,N f ) 
[2rh) N fNf\ W Nf {0,l,...,N f -l)^ 



^'Vo = h ^„,«J m 'i ' ' ,;j;^ v,iv /+ i(a,m). (C.7) 
After substituting (C.5) we find 



( e ? !0 )m-o = L im m v je ViN +l {a, rfi) ■ (C.8) 

The chiral limit of J€ vk {a, m) has to be derived differently for v = and for v ^ 0. For v = and m —> 
only the first integral in Eq. (4.34) contributes as the A-term is absent and the second integral vanishes 
because the integration range is empty. Hence, 

g— 2d noo u 2 

je ak {a,Q) = / du{-) k u k+l e-*sK k {u). (C.9) 

4a Jo 

The integral in Eq. (C.9) can be evaluated analytically in terms of an incomplete gamma function. As- 
suming that n is a nonnegative integer and Rez > 0, we find from Ref. [37, Eq. (6.631.3)] 

J n (z)= [ duu n+l e~^K n {u) = 2 n n\z'^e zlz W_n ± in{z), (CIO) 
Jo 2 ' 2 



where W\ )fl is a Whittaker function with integral representation [37, Erratum of Eq. (9.222.1)] 

z ti+\ e -zl2 

r(M-A+i) 



poo j j 

Wa u 0z) = / dte~ zt t^~ A ~2{l+ t)V +A ~2 . (C.ll) 

r(u-A+i)Jo 
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(C.12) 



Setting A = —[n + l)/2 and /i = nl2'm this equation and substituting in Eq. (C.10) gives 

J-oo f-n 
dte~ zt = 2 n z n+l n\e z Y{-n,z), 
o 1+ 1 

where the last equality follows from [37, Eq. (3.383.10)] and the incomplete gamma function is given by 

(C.13) 



T(a,z) = J 

After substituting the integral (C.12) in Eq. (C.9) we find 



dtt a - l e- 1 . 



je Qk {a,Q) = (-) k 2 k ~ 1 k\(2a) k T(-k,2&) . 



Substituting the latter in Eq. (C.8) with v = yields 

<ef 9 > v=0 ,m=o = + m2&) N f +1 Y{-N f -l,2&) . 



(C.14) 



(C.15) 



For v ^ we first analyze the two integrals in Eq. (4.34) for fh — * 0. For small argument m and v > the 
7-Bessel function goes to zero as given by Eq. (C.l), and hence the first integral in Eq. (4.34) vanishes 
because the integrand goes to zero like m v . The second integral in Eq. (4.34) trivially vanishes because 
the integration range is empty. Therefore, the limit of the average phase factor for m — ► is completely 
determined by the new A-term (4.33). It is interesting to note that the contributions to the phase of the 
determinant in the chiral limit originate from different terms for v = and v ^ 0. Equation (C.8) becomes 



-2a 



v>0,m=0 ' 



<+7sv-l (y-l-i)\(y + N f - j)l jm 2 ]' 
J™ E 



{v + Nf)l m-o ; .V (v- 1- i- jy.iljl 



8a 



(2d) 



In this double sum only the terms with i = will contribute when m — » so that 



v>0,m=0 



= e 



]= (v)jv /+ i ; ! 

where we introduced the Pochhammer symbol 

{a) n = a{a+ l)---(a+ n— 1) with {o)q — 1 



r! 



(C.16) 



(C.17) 



(C.18) 



to simplify the notation. For v = 1 we notice the intriguing fact that the chiral limit of the phase factor is 
independent of the number of flavors, i.e., (e 2ie ) v=1 lfl=0 - e~ 2a . 

Eq. (C.17) can be expressed in terms of incomplete gamma functions. To do so we note that using 
the Vandermonde convolution [38] 



{a+b) n = £ 

fc=0 



{a) n - k {b) k 



and the identity 

C-fl)» = (-l) B (a-n + l)», 
we can rewrite the coefficients of Eq. (C.17) as 



Cv-j);v /+1 N £ k < 



N f + 1 
k 



U-k+ih 



(y+N f + l-k) k 



After substituting this expression in Eq. (C.17) we find 



-2d 



v>0,m=0 



i(jvy+i,-v-u 

E 

k=0 



(-)' 



JV/+1 
k 



{v + N f + \-k) kFk 



I(j-Hl) t — 



(C.19) 



(C.20) 



(C.21) 



(C.22) 
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Because [j - k + 1)^ = for k > j only terms with k < j contribute to the sum in Eq. (C.21). Therefore 
the second sum in Eq. (C.22) starts with j - k. Moreover, because j < v — 1 only terms with A; < v - 1 
contribute, which explains the upper limit of the first sum. Using the series expansion of the incomplete 
gamma function with positive integer first argument [37, Eq. (8.352.2)] , 



n-1 z j 

T{n,z) = {n-\)\e- z Y; -tt 



Eq. (C.22) can be simplified to 



min(JVy+l,v-l) 



E 



v>0,m=0 (y + N f )\ 



(-)' 



N f + 1 
k 



(v - k) Nf+ i{2a) k Y{v - k, 2d) . 



(C.23) 



(C.24) 



Because (v - k) jvy+i = for k > v the sum in Eq. (C.24) can be extended to Nf+1 for any v > 1. In this 
case, for n < Eq. (C.23) should be replaced by the usual integral definition (C.13) of Y{n,z). Thus 



Nf+1 



v>0,m=0 ' 



(v + NfV. fc f 



E (-) fc 



N f + 1 



(v-k) Nf+ i (2d) K T (v - fc, 2d) 



(C.25) 



for any Nf and v > 1. Even though Eq. (C.25) was derived for v > 0, it can formally be continued to v = 0. 
In this case we have (- fc) jvy+i - for k = 0, . . ., Nf so that only the term with k- Nf + 1 contributes to the 
sum. We obtain 



l™< e f y ) v >o,m=o = ^/ + m2a) N f +l r(-N f -l,2&) . 



(C.26) 



This reproduces, somewhat surprisingly, the correct v = result (C.15), even though it originates from a 
different term in Eq. (4.34). 



D. Thermodynamic limit 

The thermodynamic limit is defined by a = /j, 2 F 2 V — ► oo, m — mVL — ► oo, and = rrifVL — ► oo for 
/ = 1, . . . , ivy. As in Sec. 5.2, we assume for simplicity that m = mi = ... = m/vy • In the following we show 
that for 2alm< 1, the phase factor in the thermodynamic limit is given by 



i,\N f +l 



2a\ n f 



m 



(D.l) 



The starting point is Eq. (3.30) for the phase factor in the unquenched case with Nf dynamical quarks of 
equal masses. We now compute the thermodynamic limit of the Wronskian W n {k\,. . . , k n ). Starting from 
the definition I v ^{m) - fh Iv+fciffi) it is easy to show that the p-t\i derivative is given by 



p 

q=0 



m^'/^lm), 

{k-p+q)\ v+k 



(D.2) 



where the latter expression only contains derivatives of the 7-Bessel function. After substituting this 
expansion in the Wronskian determinant and using the asymptotic expansion 



I v (m) = 



i+E 



y=i mi 



+ e- m {-) 



(D.3) 



of the 7-Bessel function one can show, using basic properties of determinants, that in the thermodynamic 
limit only the £7 = term of Eq. (D.2) contributes to leading order so that 



W n {k l ,...,k n ) 



I2n 



m^ ki - n2 ' 2 A n {h,...,k n ) 



(D.4) 
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for arbitrary {A;, }, where A n (jfci, . . . , k n ) is a Vandermonde determinant. Using this expression and Eq. (C.6), 
the Wronskian in the denominator of Eq. (3.30) becomes 



Nr Nf-1 



(D.5) 



In the thermodynamic limit J€^ k {a, m) can be computed using Eq. (C.2) of Ref. [28], resulting in 



V 2m \ m , 



(D.6) 



which is independent of v. To compute the thermodynamic limit of the Wronskians in Eq. (3.31) we use 



JV f +l 



A Nf+1 (0,...,k-l,k+l,...,N f + 1) = 
Nf+1 * (Nf + l-ky.kl 



and £ f= -(Af / + 2)(AT / +1)-A; (D.7) 

!=0 2 



to show that 



W N +i(0,...,k-l,k+l,...,N f + l) 



\ N f + 1 N f + 1 

m 2 



: il ^'=l £ - 

(jy^+i-fc)!fc! 



After substituting Eqs. (D.6) and (D.8), the thermodynamic limit of Eq. (3.31) becomes 



1 / e n 



m \ N f N f 



(Nf+1 \ 



W Nf (a, m) ~ - 1 -= I m * 
From the binomial theorem we know that 



Nf+1 



1 



1 



V /=i J to Wf + l-kW\ m 



(1 + x)" 



4a) k 



£ {n-k)\k\ n\ 



and hence 



W Nf {a,m)~2 N f 



e m \ N f N f 



In) 



m 2 



(N f \ 



*,\N f +l 



1 



2a V v / 



(D.8) 



(D.9) 



(D.10) 



(D.ll) 



After substituting Eqs. (D. 1 1) and (D.5) in Eq. (3.30) we find the thermodynamic limit for the phase factor, 

\N f ~/ ( N f „,\ ,/>^Nf+l 



(e?T = 



1 ^/(^^(ii^jd-f)^ 



WW (n^HA) 

which simplifies to (1 - 2alm) N f +1 as given in Eq. (D.l). 



N f 



(D.12) 



E. Numerical random matrix simulations 



In this appendix we describe the numerical simulations of random matrices used to verify the analytical 
results derived in the main body of the paper, for both trivial and nontrivial topology. This procedure 
also illustrates the potential usefulness of numerical simulations in cases where analytical results would 
not be immediately accessible. 

We performed numerical simulations of random matrices in the chiral GUE with chemical potential 
for the quenched case. As mentioned in Eq. (2.1), these random matrices can be constructed as 



DQi) = 



/<P + J u v T 
/<P t + ^T' + 



(E.l) 
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where $ and *P are complex [N+v) x N matrices generated according to the Gaussian weight function 
w{X) oc exp(-NtrX + X) = exp ( - N^ \X ke \ 2 ) = U ex P ( " W(ReX fcf ) 2 ) exp ( - N(\mX ke ) 2 ) . (E.2) 



kC 



kC 



The last expression shows that the real and imaginary parts of each matrix element are i.i.d. random 
numbers drawn from the Gaussian distribution 



w{x) oc exp (-Nx 2 ) 



(E.3) 



with standard deviation 1 / V2N. 

As we want to investigate the microscopic limit of the theory, we will keep a = 2Na and m = 2Nm 
constant, while taking N large enough to approach the microscopic limit, in which N — » oo. Hence, 
when generating and diagonalizing the matrices from Eq. (E.l), the chemical potential will be scaled as 
[i - \J &/2N. Each random matrix Difi) is then diagonalized and the real part of the phase factor of its 
determinant computed with 



cos20 = cos 



2JV+v 

2 V arg(A,- + m) 
l=i 



(E.4) 



where the A, are the eigenvalues and m = m/2N. For a sample with N s random matrices the real part of 
the average phase factor will be given by 



cos20(v, a, m) — — V cos20; , 



(E.5) 



where 6j is the phase of the determinant of the y'-th random matrix in the sample, given by Eq. (E.4). For 
simplicity we have omitted the subscript 5 (for the microscopic limit) on cos 28. The average of the imag- 
inary part will be zero within the statistical error because of the symmetry properties of the ensemble 
and is therefore disregarded in our analysis. 

The chiral symmetry of the matrix can be used to improve the efficiency of the computer implemen- 
tation by transforming the (2JV+ v) x (2N+ v) diagonalization problem to one of size N x N. Let us first 
rewrite Eq. (E.l) as 

'0 A 
B 



(E.6) 



The eigenvalue equation 



can then be written as 



Di/fiv- Av 



fo a) 












[b 0j 






= A 




■-{ 




v 2j 




v 2j 





Av2 = Xvi , 
Bv\- Xv2 



(E.7) 



(E.8) 



with complex eigenvalues A and eigenvector decomposition v — {v\,V2), where V\ and v 2 are complex 
vectors with [N + v) and N elements, respectively. Without fine-tuning D{[i) has exactly v zero modes 
that obey the two homogeneous linear systems 



Av 2 = , 
Bvi = 0. 



(E.9) 



The first system contains N + v linear equations with N variables, and the second one N equations with 
N+v variables. If both A and B are of rank N (no fine-tuning), the first homogeneous system only has 
solutions v 2 - 0, while the second system is underconstrained and has v linearly independent solutions 
for v\. Hence the zero modes will be represented by v eigenvectors (fifc,0) with k — l,...,v. 
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Moreover, it is clear from Eq. (E.8) that each nonzero eigenvalue A with eigenvector (v\, vz) will be 
paired with an eigenvalue - A with eigenvector {v\, -i>2l- This is a consequence of the chiral symmetry of 
the problem. 

These properties of the spectrum of D{ji) will now be used to transform the diagonalization problem 
from order 2N+ v to order N. For this, let us first multiply Eq. (E.8) from the left with D{/f), 



D z in)v = 



AB 
BA 



yV 2 j 



(E.10) 



where AB is an {N + v) x [N + v) matrix and BA has dimension Nx N. Clearly, each nonzero eigenvalue 
A 2 of AB (with eigenvector V\) is also an eigenvalue of BA (with eigenvector v-i). However, AB has v 
additional eigenvalues, which necessarily correspond to the zero modes (^1^,0) satisfying Eq. (E.9). 

This can be used to expedite the numerical simulations. It suffices to diagonalize the N x N matrix 
BA to find the N nonzero eigenvalues A 2 . We then know that the eigenvalues of D(/i) are the N pairs 
(Ajj-Aj) supplemented by v eigenvalues equal to zero. The determinant for a fermion of mass m will 
then be given by 



N 



det(DQu) + m) = m> '\\{m 2 - A?) , 

2=1 



and the real part of its phase factor is 

cos 28 = cos 



2£arg(m 2 -A 2 ) 



(E.ll) 



(E.12) 



which replaces Eq. (E.4) in our practical simulations. Note that the cost of the additional multiplication 
of an AT x (N+ v) by an (AT + v) x N matrix to construct the product B A is negligible compared to the cost 
of the diagonalization. 
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